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Michael Reichert and Holger Stark 

Fachbereich Physik, Universitat Konstanz, D-78457 Konstanz, Germany 
E-mail: michael . reichert@uni-konstanz . de 

Abstract. Particles suspended in a viscous fluid circle in optical vortices generated by 
holographic optical-tweezer techniques [Curtis J E and Grier D G 2003 Phys. Rev. Lett. 
90 133901]. We model this system and show that hydrodynamic interactions between the 
circling particles determine their collective motion. We perform a linear-stability analysis to 
investigate the stability of regular particle clusters and illustrate the limit cycle to which the 
unstable modes converge. We clarify that drafting of particle doublets is essential for the 
understanding of the limit cycle. 

PACS numbers: 05.45.Xt, 47.15.Gf, 47.85.Np, 82.70.Dd, 82.60.Yz 

1. Introduction 

Hydrodynamic interactions occur between particles or bodies whenever they move relative 
to each other in a viscous fluid. Due to their long-range nature, they are important for the 
dynamic properties of colloidal suspensions exemplified by the self- and collective diffusion 
Q]|2]|3]|4), sedimentation |3]|6|[7), or the aggregation of particles 1 8 1. Through the correlated 
motion of a pair of colloids trapped in optical tweezers, one can directly measure the effect of 
hydrodynamic interactions ll9l ll0lll 1111 21 . Furthermore, hydrodynamic interactions give rise 
to interesting collective behavior, e.g., periodic or almost periodic motions in time 1131 1 141 
or even transient chaotic dynamics in the sedimentation of particles [15|. And they lead to 
pattern formation of rotating motors 1161 with a possible 2D melting transition of biological 
motors like ATP-synthase embedded in a membrane 1171 . Hydrodynamic interactions are 
treated in the low-Reynolds-number regime which is also relevant for biology. It determines 
the problem of how microorganisms move forward [18]. Certain bacteria accomplish this, 
e.g., by cranking helical flexible rods 1 19 1. In addition, recent experiments show that laminar 
flows initiate the asymmetries between the right- and left-hand side of the body at an early 
stage of the embryonic development I20l l21ll22"l . 

The work presented here is motivated by experiments of Curtis and Grier |23|. They 
created toroidal optical traps known as optical vortices with the help of holographic 
techniques. In the bright circumference of an optical vortex, they could trap particles 
which circulated around the ring due to scattering forces which result from the orbital 
angular momentum of light. In this article, we model the system of circling particles. We 
demonstrate that hydrodynamic interactions determine their interesting collective motion 
which we analyze by methods used in the study of non-linear dynamics. In particular, we 
identify a limit cycle which is governed by drafting of particle doublets. 

The outline of the article is as follows. In section [2] we first define our model system 
of particles moving in a ring-like trap and summarize the basic equations for the description 
of the dynamics. The stability of regular A^-particle clusters in the ring is investigated in 
section |5] by means of a linear stability analysis. Section |4] then describes the nonlinear 
dynamics of perturbed clusters. We introduce the periodic limit cycle of the three-particle 
system and perform a harmonic analysis. The influence of the trap strength on the dynamics 
is discussed. Furthermore, we briefly address the dynamics in fairly weak traps. Finally, in 
section|5] we present the particle velocities as a function of ring radius and particle number. 
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Figure 1. Definition of our model system of an optical vortex. The particles are driven by a 
constant tangential force and trapped in the ring by a harmonic radial force. In order to avoid 
overlaps, the particles repel each other at very short distances. 



2. Modelling particle dynamics in a toroidal optical trap 

We consider non-Brownian, equal-sized spherical particles suspended in a viscous fluid in 
the regime of low Reynolds numbers whose mutual interactions are solely of hydrodynamic 
origin. We mimick the basic features of particles captured in an optical vortex by applying a 
constant tangential force to each particle and by keeping the particles on a ring of radius 
R by means of a harmonic radial trap with force constant K r (figure^. The particle motion 
is effectively two-dimensional in the plane of the ring (z = 0); thus, the particle positions 
are best described by polar coordinates, and fa. With the radial and tangential unit vectors 
at the position of particle i, e\ = (cos fa, sin fa, 0) and ef = (— sin fa, cos fa, 0), the total 
external force acting on particle i then reads 



Note that the assumption of a constant tangential driving force is a rough 
approximation since the intensity profile of an optical vortex is modulated along the ring 
as described by the topological charge £ of the vortex 1231 . However, if i is large so that 
the period of the modulations is small compared to the particle size, then our approximation 
seems to be reasonable. 

In order to prevent that particles overlap in the simulations, we add a repulsive interaction 
which becomes relevant only when two particles are very close to each other. We use a hard- 
sphere-type interaction potential where a "hard core" with diameter 2a is surrounded by a 
"soft" repulsive potential (figure^ |28 1: 



Here, is the center-to-center distance between particles i and j, and a is the particle radius. 
This choice is primarily due to numerical reasons, but can nevertheless be justified from the 
physical point of view since even hard-sphere-like colloids show a "soft" repulsive potential 
at very short distances 1281 . We choose the commonly used Lennard- Jones exponent m = 12 
so that the repulsion only acts at very short distances, and we tune the prefactor W such that 
the minimal inter-particle gap — 2a in the simulations is of the order of 10~ 4 a. 

In the regime of low Reynolds numbers, the flow of an incompressible fluid with 
viscosity 77 obeys the Stokes or creeping flow equations 1241 1251 l2l 



Fi = F^ef — K r (r.i — R)e\ . 



(1) 




(2) 



ryV 2 w ~X7p = 0, V • u = , 



(3) 
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where u is the flow field and p the hydrodynamic pressure. Imposing stick boundary 
conditions on the surfaces of all particles suspended in the fluid, the motions of the particles 
are mutually coupled via the flow field. Due to the linearity of 0, the translational and 
rotational velocities of the particles, V i and Wj (i = 1, . . . , AT), depend linearly on all external 
forces and torques acting on the particles, Fj and Tj [24 251121: 

N N 

^ (//I'/-',- /<;r ; ) . m ^(//;;/<\ • ,r7\) . (4) 

i=i j=i 

The central quantities are the 3 x 3 mobility tensors fi^, fif-, and /x-j. They depend on 
the current spatial configuration of all particles, i.e., the set of position vectors {r\, . . . , rjy} 
in the case of spherical shape. 

As there are no external torques in our model, i.e., Tj = 0, and as we are not interested 
in the rotational motion of the particles, the remaining equation of motion for our problem is 

N 

ri = v i = J2t4 t 3 F J- ( 5 ) 

3=1 

Since the mobility tensors /x** are nonlinear functions of all particle positions {r^}, 
describes the coupled nonlinear dynamics of N particles. 

The first-order approximation (with respect to inverse particle distances ry) for the 
mobilities of particles moving in an unbounded and otherwise quiescent fluid is the well- 
known Oseen tensor 1 2 1 

where /i* = (67rr/a) _1 is the Stokes mobility for a translating sphere; the self-mobilities are 
fAl = A**!- All the other mobility tensors in (0} vanish in this approximation. Note that 
the Oseen tensor is the Green function of the Stokes equations l|3}, i.e., it considers pointlike 
particles and hence does not include rotational couplings. 

Higher order approximations of all mobility tensors can be calculated, e.g., via the 
multipole expansion method |26|. It has been implemented in the numerical library 
HYDROLIB [27 1, which we use in our simulations. For the numerical integration of the highly 
nonlinear equation of motion l|5}, we apply a fourth-order Runge-Kutta scheme. The time 
step is chosen such that the corresponding change in the angular position is of the order of 1°. 



3. Stability analysis of regular clusters 

In this section, we study the stability of regular iV-particle clusters. By "regular" we mean 
configurations with iV-fold symmetry, like the ones shown in figure|2] Due to this symmetry, 
all particle positions are equivalent, i.e., the angular velocity <f>i has to be the same for all 
particles. The radial positions of the particles do not change either, which can also be justified 
by a simple symmetry argument. Let us start with a regular configuration where = R for 
all particles, so there is not radial force component in Q. Due to linearity, oc F^, which 
means that changing the direction of F^ inverts the direction of f^. However, since there is 
no difference if the cluster rotates clockwise or counterclockwise, the radial velocity must be 
the same for both directions, i.e., = and = R. Therefore, a regular cluster remains 
unchanged and rotates with constant frequency fijy. 

In regular clusters, the particles are well separated for sufficiently large radius R, and 
we may approximate the hydrodynamic interactions by the Oseen tensor (jfji. The particle 
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Figure 2. Regular TV-particle clusters. 



positions of the regular clusters are given by = R and fa = = (2n/N)i + flj^t. With 
Fi = F&ef, we derive from (|5j the rotational frequency of an iV-particle cluster, 



1 



3^1 1 



3cos$ 



n x = 



R 



(7) 



where = \/2{R/a)^/\ — cos$ i:) and $.y = (27r/iV)(j — i) are the respective spatial and 
angular distances between particles i and j. 

To investigate the stability of these iV-particle clusters against small radial and angular 
particle displacements, we introduce T{ = R[l + Spi(t)] and fa — $i(i) + Sfa(t). Using the 
forces {0 and the Oseen tensor (|6j, we linearize the equation of motion (|5} in terms of the 
small perturbations Spi and 6 "fa: 
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where r = ilii is the reduced time, and the dimensionless parameter K. = K r RjF^ 
measures the strength of the radial trap relative to the tangential driving force. The 
corresponding eigenvalue problem was solved numerically by using the computer-algebra 
package MAPLE. This analysis reveals that there are the following types of stable and 
unstable eigenmodes (depending on the number of particles) for the coupled displacements 
(Spi(t), . . . , 6pn(t), <50i(t), . . . , S<Pn(t)): (i) a constant angular shift, (ii/iii) non- 
oscillating damped or unstable, and (iv/v) oscillating damped or unstable. Figure [3] shows 
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Figure 3. Representative example of the angular coordinate of an unstable oscillating mode 
(here: N = 4, Ft/ a = 10, ?C = 10). The linear theory and the simulations agree very well up 
to amplitudes of several percent of 2n/N. (Only every 100th data point from the simulation 
is plotted.) 



(i) (ii) 

N 

A = A e R_ 


(iii) 

A G M+ 


(iv) 
Re A < 


(v) 
Re A > 


odd 1 1 

even 1 2 


1 


AT- 1 
2V-2 


N - 1 
N - 2 



Table 1. Number of eigenmodes classified by their eigenvalues A for even and odd particle 
numbers N. M± means a positive or negative real number, (i) constant angular shift, (ii/iii) 
non-oscillating damped or unstable modes, (iv/v) oscillating damped or unstable modes. 



an example of an unstable oscillating mode (type v). In table [2 we give the number of 
eigenmodes classified by their eigenvalues A for even and odd particle numbers. 

The occurence of stable and unstable modes can be compared to a saddle point in the 
framework of the analysis of dynamic systems. A cluster configuration with arbitrary small 
radial and angular displacements of the particles is practically unstable since contributions of 
unstable modes will grow, while stable modes will relax to zero. 

For N = 3, we have also determined the eigenvalues analytically by an expansion to first 
order in a/R: 

Ai =0 (type i) , 

y/3 a 

To R, 



A, 



-K[ 1 



(type ii) . 
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(type iv) , 



(type v) . 



In this case, a non-oscillating unstable mode (type iii) does not exist; it only occurs for even 
N (see table [TJ. 
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Figure 4. Dynamic transition in the angular coordinate from an unstable linear mode to the 
periodic limit cycle for three particles (R/a = 10). The main curve and the insets (a), (b), and 
(d) show the angle (j>(t) relative to fijvt as a function of time. In the inset (c), <j>(t) — (Sl)t 
is plotted. Plots (a)-(c) correspond to a stiffness of K. = 10. For comparison with (b), the 
transition regime for K = 100 is shown in (d). The dashed lines in (b) and (d) indicate the 
limit cycles. 
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Figure 5. Basic mechanism of the periodic limit cycle. Two particles in close contact move 
faster than a single particle. When this pair reaches the third particle, they form an intermediate 
triplett, and finally, the first two particles escape as a new pair. 



4. Nonlinear dynamics of perturbed clusters 

The amplitude of an unstable mode grows up to a certain magnitude in agreement with the 
linear analysis, until the nonlinear dynamics takes over (see figure^. The amplitudes saturate, 
and the system finally tends towards a periodic limit cycle with oscillating particle distances. 
In figure|4]we show an example of such a dynamic transition from an unstable linear mode to 
the periodic limit cycle by plotting the angular coordinate cf)(t) relative to fi^i as a function 
of time. Note that il^ is the rotational frequency of the regular N -particle cluster introduced 
in0. 

The mean slope in each of the two regimes (linear mode and limit cycle) gives a well- 
defined mean orbital frequency (f2), which is Sl^v for the linear mode. Clearly, (fl) and 
therefore the mean velocity of the particles is increased by the transition to the limit cycle, 
which means that the mean drag force on the particles is reduced. 

The character of the transition depends on the trap stiffness. It is quite sharp for weaker 
traps, where the transition takes place within a few oscillations, as shown in figureUb). For 
stronger traps, the mean frequency increases smoothly from the linear regime to the limit 
cycle, which is illustrated in figure|4jd). Furthermore, the onset of the transition is shifted to 
later times when the trap becomes stiffer. 
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Figure 6. Fourier spectrum [cos -coefficients a(ui) and sin-coefficients b(u>)] of the periodic 
limit cycle of three particles. Note that the frequencies ui are defined relative to the mean 
orbital frequency (SI). 
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Figure 7. Characteristic frequency u>* of the periodic limit cycle of three particles versus ring 
radius R. The different curves are related to different trap stiffnesses K r . 



Figure [5] introduces the basic mechanism underlying the periodic limit cycle. Two 
particles in close contact move faster than a single particle since the friction per particle is 
reduced due to the well-known effect of drafting. When such a pair reaches the third particle, 
they form a triplett for a short time. In this configuration, the mobility of the middle particle is 
larger since it is "screened" from the fluid flow by the outer particles. It therefore pushes the 
particle in front, so that the first two particles "escape" from the last one. The same principle 
also holds for more than three particles. In these cases, there can be more than one pair of 
particle. 

For the limit cycle of three particles, we have performed a harmonic analysis using 
Fast Fourier Transformation. In order to separate the fast orbital dynamics of the particles 
(characterized by the mean angular frequency (51)) from their relative motions, we calculate 
the Fourier transform of <p(t) — (f2)i = ^ w [a(u;) coswi + b(w) sinwt]. This yields the 
fingerprint of the dynamics relative to the mean circling velocity. In figure [6] we show the 
corresponding Fourier spectrum. Besides the characteristic frequency oj* of the limit cycle, 
the higher harmonics 2ui* and 3uj* are also very pronounced. 

The characteristic frequency lj* decreases with increasing ring radius, as shown in 
figure[7] because it takes longer for a particle pair to reach the third particle on a larger ring. 
At constant ring radius, the characteristic frequency increases with the trap stiffness since 
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Figure 8. Dynamics of four particles in a weak radial trap. The motion shown is relative to 
the mean rotation along the ring trap. Only the relevant part of the ring is shown, the center 
of the ring is located at the bottom of the snapshots. The particles are running from right to 
left. During one cycle, the first and the third particle are exchanged. The particles are shaded 
differently in order to distinguish them. 



particles in a stronger trap are better aligned along the ring, which makes the mechanism of 
the limit cycle more effective. 

If the strength of the radial trap is decreased, the radial displacements of the particles 
increase so that they can pass each other. This happens at a reduced trap stiffness JC = 
K r R/F^ of the order of 1 . For three particles and four particles, the limit cycle then consists 
of a compact triangular or rhombic-shaped cluster circling and rotating in the trap. 

For four particles, we also find the limit cycle illustrated in figure[8] It occurs if the radial 
trap has medium strength so that compact particle clusters only have a finite life time. The 
first particle on the left is pushed in outward radial direction by the particles behind it and 
subsequently passed by the second particle. It, then, pushes the third one in inward radial 
direction, which in turn takes over the lead of the chain. Finally, the initial state is rebuild, 
where the first and the third particle are exchanged. Note that there is an intermediate state 
where the particles form a compact rhombic-shaped cluster, which, however, is not stable. 

5. Particle velocities 

In figure[9](left), we study the circling frequencies of regular clusters as a function of particle 
number N for different ring radii R. At constant radius, the frequency increases with the 
number of particles since then the particles are closer together which reduces the drag. Note 
that for less than four particles, the frequency is reduced relative to the single particle value. 
Here, hydrodynamic interactions across the ring increase the drag on the particles. This effect 
is clearly more pronounced for small rings. 

At constant particle distance 2ttR/N or constant line density of the particles N / (2ttR), 
the drag is reduced when the ring radius increases [see figure[9](right)]. This is due to the fact 
that in a ring with smaller curvature the particles are better "aligned" behind each other. 

In the limit cycle, discussed in the previous section, the mean orbital frequency is larger 
compared to the velocity of the corresponding regular configuration. The quantitative effect is 
illustrated in figure ^3 We attribute this to the drafting of particle doublets which obviously 
reduces their drag force. We observe that the dependence of the limit-cycle velocity on the 
trap stiffness is weak, in contrast to the characteristic frequency of the angular displacements, 
as already discussed in figure^ 

6. Conclusion 

We modelled the circling of particles in an optical vortex and illustrated how hydrodynamic 
interactions govern the nonlinear dynamics of the coupled particle motion. We hope that our 
theoretical investigation initiates a detailed study within experiments. Possible extensions of 




Figure 9. Circling frequencies of regular ^-particle clusters, normalized to the frequency of 
a single sphere, as obtained from simulations. The two plots show the same data, left: as a 
function of particle number, right: as a function of particle line density. The arrow gives the 
direction of increasing ring radius. 
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Figure 10. Mean circling frequencies of three particles in the periodic limit cycle, normalized 
to the frequency of a single sphere. The upper three curves are plotted for different trap 
stiffnesses K r . For comparison, the lower curve shows the values for the regular three-particle 
cluster. The errorbars indicate the range within which the orbital velocities vary in the limit 
cycle. 



our work concern the tangential driving force which could be modulated along the ring or 
which could possess stochastic contributions. 
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